On the critical behavior of a lattice prey-predator model 
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The critical properties of a simple prey-predator model are revisited. For some values of the 
control parameters, the model exhibits a line of directed percolation like transitions to a single 
absorbing state. For other values of the control parameters one finds a second line of continuous 
transitions toward infinite number of absorbing states, and the corresponding steady-state exponents 
are mean-field like. The critical behavior of the special point T (bicritical point), where the two 
transition lines meet, belongs to a different universality class. The use of dynamical Monte-Carlo 
method shows that a particular strategy for preparing the initial state should be devised to correctly 
describe the physics of the system near the second transition line. Relationships with a forest fire 
model with immunization are also discussed. 
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I. INTRODUCTION 

The study of prey-predator systems has attracted a lot 
of attention since the pioneering works of Lotka [Q and 
Volterra Q . Working at a mean- field level (homogeneous 
populations) they showed that, depending on the initial 
state, the system can evolve toward a simple steady-state 
or a limit cycle, in which the populations oscillate peri- 
odically in time. 

An important question is the understanding of the role 
played by the local environment on the dynamics (spa- 
tial effects) H and, accordingly, many extended prey- 
predator models have been studied during the past years 
|^-^| . Recently, a simple prey-predator model was intro- 
duced by some of us |10(| . Although governed by only 
two control parameters, this model exhibits a rich phase 
diagram. As a function of the two control parameters A a 
and Ab, which are the growth rates of prey and preda- 
tor respectively, two different phases are observed: a 
pure prey phase (P), and a coexistence phase of prey 
and predator in which an oscillatory (O) region and a 
non-oscillatory (NO) region can be distinguished. For a 
system size L — > oo, these three different domains meet 
at a particular point, called T = (X^,XT) (precise defi- 
nitions are given below). It was shown |Hj that = 
and Xf ks 5.0 ± 0.3. For A a > 0, a phase transition 
line between the pure prey phase and the coexistence 
phase is present, and the critical exponents along this 
line are the ones of directed percolation (DP) [ pT[ . How- 
ever, it was also observed that when the growth rate of 
prey A a — > and A;, > X[ , the model undergoes a non- 
DP continuous phase transition. Since the directed per- 
colation is a generic universality class for models with 
absorbing states (unless some special conditions are sat- 



isfied |T^1), existence of such a transition is certainly 
surprising. These two lines of different continuous non- 
equilibrium phase transition meet at the bicritical point 
T Jl3| , and one forecast that the critical behavior at this 
particular point may also belong to a new universality 
class. 

The goal of this paper is to study in more details the 
properties of these non-DP phase transitions. First, we 
perform extensive steady-state simulations, which con- 
firm the non-DP character of the transition in the limit 
A Q — » and A& > A^. 

As it was already shown ||lo|| , in this limit the model 
exhibits oscillatory behavior. But, in addition to that, for 
A Q = the model has infinitely many absorbing states. 
These two properties are responsible for a rather peculiar 
behavior of the model, which becomes particularly trans- 
parent when the model is examined using the dynamical 
Monte Carlo method. When applied to models with in- 
finitely many absorbing states, this dynamical method 
uses the so-called natural absorbing states, which are the 
most likely states to be reached by the dynamical evolu- 
tion of the system. We show, however, that this common 
but somehow heuristic procedure fails here. Indeed, for 
the present model, natural absorbing states contain only 
short-ranged islands of prey on which spreading is not 
critical. To restore criticality of spreading we generated 
the absorbing states using a quasi-static approach to the 
critical point. This example shows that for some models 
with infinitely many absorbing states a special approach 
is needed to examine the dynamical properties of the crit- 
ical point. 

The paper is organized as follows. In Sec. |l[ the 
model is defined and some of its properties discussed. A 
thorough investigation of the critical behavior has been 
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made using two different complementary approaches. In 
Sec. I_0, the critical behavior is investigated using steady 
state properties while in Sec IV one uses dynamical 
Monte Carlo method. 



and Ah - 



It is shown that for A a — > 
the steady-state exponents are indeed 



mean-field-like, while the dynamical exponents are non- 
universal depending continously upon Af,. Nevertheless, 
trace of the mean-field character of the transition shows 
up in scaling relation among dynamical exponents. The 
critical behavior at point T is also investigated, and it 
turns out that the corresponding exponents belong to a 
new universality class. Finally, physical arguments ex- 
plaining the above findings are given in Sec. |Vj 



II. MODEL 

The model used in Ref. is defined as follows. Each 
cell of a two-dimensional square lattice (of size L x L, 
with periodic boundary condition), labeled by the in- 
dex i, can be at time t, in one of the three following 
states: at = 0, 1, and 2. A cell in state 0, 1, or 2 corre- 
sponds, respectively, to a cell which is empty, occupied 
by prey, or simultaneously occupied by prey and preda- 
tors. The transition rates for site i are (i) — > 1 at rate 
A (ni,i +ni >2 )/4, (ii)l — > 2 at rate A b (n ii2 )/4, and (iii) 
2 — > at rate 1, where n i tT denotes the number of near- 
est neighbor sites of i which are in the state a. The first 
two processes model the spreading of prey and predators. 
The third process represents the local depopulation of a 
cell due to overly greedy predators. The rate of the third 
process is chosen to be 1, which sets the time scale, hence 
t, as well as A a and A& are dimensionless quantities. 

The properties of this model have been investigated 
both by mean- field and Monte-Carlo methods |icfl . The 
Monte-Carlo result extrapolated to the case L — > oo is 
summarized in Fig. |l|. The transition line between the 
prey phase and the coexistence phase, Ai*(A a ), belongs 
to the directed percolation universality class JpTj, as ex- 
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pected, and terminates at the point T = (A„ 
where the P, O, and NO domains meet. For A;, > X[ , the 
transition between the oscillatory domain of the coexis- 
tence phase and the prey phase takes place at A Q = 0. 
Along this transition line the predator density approaches 
zero as a power law b ~ Af 2 , with fa ~ 1, and so, does not 
belong to the DP class, which is somehow unexpected. 
The value fa ~ 1 lead to the conjecture JlQ] that this 
second transition could be mean-field-like. There is a 
crossover between the O and NO parts of the coexistence 
phase. The purpose of the present study is to give a com- 
plete description of the nature of the transitions near the 
line A a = for A& > . 

It is worth to mention, that our model is closely related 
to a model introduced by Drossel et al. |m| to investigate 
the effect of immunization in an extension of the simple 
forest- fire model [ pT5[ . This three-state model (0: empty 
site, 1: tree, and 2: burning tree) differs from our model 



in some details: the growth rate of a tree (a : — > 1) 
is p, independently of the environment, and a tree is ig- 
nited (a : 1 — ► 2) with rate (1 — 3)0(712), (© is the 
usual Heaviside- function) . This second process models 
the immunization of trees against fire. The third process 
(<7 : 2 — > 0) happens at rate 1. For non zero immunity 
and p > Albano showed that the transition toward 
a single absorbing state is DP like, while for p = (at 
the end point of the DP transition line), the transition 
belongs to the dynamical percolation universality class, 
and the absorbing state is not unique. 
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FIG. 1. Phase diagram obtained by extrapolation of the 
simulation results to the case L — > 00. The solid line repre- 
sents the DP transition between the prey phase (P) and the 
non-oscillatory part (NO) of the coexistence phase ( The • 
symbols are the simulated values while the solid line is just 
a guide to the eyes). The □ symbols delimit the crossover 
between the oscillatory (O) and non-oscillatory (NO) regimes 
present in the coexistence phase. The arrows correspond to 
the path described in the text, along which the critical expo- 
nents have been measured. 



III. STEADY STATE STUDY OF THE CRITICAL 
BEHAVIOR 

Extensive Monte-Carlo simulations for system sizes up 
to 4000 x 4000 have been performed to investigate the 
behavior of the predator density b, for A a — > and three 
different values of Ah, namely, A& = 4.67, 5.0, and 6.0, 
following trajectories of type 2 and 4 in Fig. |l|. The 
value A a = 0, A& = 4.67 corresponds to the best deter- 
mination of the end point T, obtained by the dynami- 
cal approach described below. Owing to the oscillatory 
behavior near the critical line (A a = 0, A/, > the 
system very easily evolves into an absorbing state, where 
the predators are extinct, therefore careful initialization 
is needed in the simulations. Usually 10 4 MCS "thermal- 
ization" were applied following 10 4 MCS "approaching" 
time where A a was decreased continuously. The densities 
and the fluctuations of prey and predators were averaged 
over ~ 2 x 10 5 MCS for each A a , Xb points. It is found 
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that b ~ Af 2 for A a — > 0. In order to see corrections to 
scaling, we compute the effective exponent 

n n lnft(A„(i))-ln&(A a (i-l)) m 

lriA a («) - mX a {i - 1) 

where A a (?) and A a (i — 1) are two consecutive values of 
the control parameter A a . 



in the active state. One considers an ensemble of trials 
starting from the same initial state. Certain dynamical 
quantities exhibit a power law behavior when the system 
is critical. For example, the survival probability behaves 
as 



P(t) ~ r 



(6) 
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FIG. 2. Predator density critical exponent 02 obtained for 
several values of A a and A;, = 4.67 (o), 5.0 (□) and 6.0 (O). 

As the Fig. || shows, for A^ = 5 and 6, the linear ex- 
trapolation of /3 e f f converges to ~ 1 within statistical 
errors 



fa(X b = 5.0) = 1.01(1), 
fa(X b = 6.0) = 0.96(2). 

For Ab = 4.67 (trajectory of type 2 in Fig. 
somewhat higher value, and one finds 

fa(X b = 4.67) = 1.33(4). 

The measurement of the fluctuations 

Xb = L 2 <(6-<6)) 2 ) ~A-" 



(2) 
(3) 

it goes to a 
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are less precise and we estimate 7 = —0.6(16) at X b = 
4.67, and 7 - for A b = 5, 6. 

For Ah > Xj these values are consistent with the pre- 
vious prediction [|l(| fa ~ 1. However, at the bicritical 
point T the value of fa is completely different and thus 
belongs to a new universality class. 



IV. DYNAMICAL STUDY OF THE CRITICAL 
BEHAVIOR 

A very useful technique to study the critical proper- 
ties of a system with absorbing states is the so-called 
dynamical Monte Carlo method [0. In this approach, 
the system is prepared in an initial state, which is one 
of the absorbing states up to one site, which is set to be 



The deviation from this power law behavior, when the 
system is off-critical, provides a very precise way to lo- 
cate the critical point. 

The number of active sites N(t) behaves as 



N(t) - f 



(7) 



while, for the mean square spreading from the origin 



R 2 {t)~t z 



(8) 



where the dynamical exponent z = 2v_l/v\\ is the ratio of 
the critical exponents of spatial (v±) and temporal (i^|) 
correlation lengths. Some scaling relations between these 
exponents can be also derived |l§| ]. 

We made various simulations using the dynamical 
Monte Carlo method and our results are summarized be- 
low. 



A. Case of A a = 



Oh 




FIG. 3. The survival probability Pit) as a function of t ob- 
tained for A a = and (from the top) Ab = 4.6, 4.65, 4.67, 
4.7, and 4.75 (trajectory of type 3 in Fig. [l]). We used system 
size L = 3000, and up to 10 5 independent runs were made for 
each value of A;,. The dotted line has the slope corresponding 
to S = 0.092. 
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First, we simulated the model on the A a = line, tak- 
ing, as an absorbing state, a lattice filled with prey (tra- 
jectory of type 3 in Fig. [[]). Measuring the survival prob- 
ability Pit), we found that A& = \J « 4.67 is the critical 
point, which separates the absorbing phase (Af, < A^) 
and the phase with annular growth (A& > Xf). Mea- 
suring the slope at Af, = Xj (see Fig. ||) we estimate 
5 w 0.095(5), which is very close to the value obtained 
for dynamical percolation, for which, in two dimensions, 
S = 0.092 Moreover, using (7) and (8) we obtained 
1] = 0.60(5) and z = 1.72(4), which are also very close to 
the dynamical percolation values. 

Note that the usual /3 exponent for dynamical percola- 
tion, which takes the value (3 ~ 0.14, is defined through 
b ~ (A? — Xb) 13 , which differs from our definition of /?2 in 



Sec. ED. Thus it is not surprising that these two expo- 
nents differ. Note also that the dynamical estimation of 
the critical endpoint Af, = A^ has been used in the static 



approach of Sec. ED 



B. Case of A a > 

The same scheme was used for A a > 0. The critical 
point was located for A a = 0.5 and 1 (trajectory of type 
1 in Fi g, pi ), and they agree with the steady state results 
of Ref. |l5J . Measuring the slope at criticality (see Fig. |i|) 
we estimate 5 « 0.45, i.e., the value compatible with DP 
p8| . The fact, that for A a > the phase transition be- 
longs to the DP universality class, was already confirmed 
using the static calculations |10|. 
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FIG. 4. The survival probability P(t) as a function of t ob- 
tained for A a = 0, Xt = 4.67 (trajectory of type 3 in Fig. |l|), 
for A a = 0.5, X b = 3.175 and for X a = 1, X b = 2.451 (tra- 
jectory of type 1 in Fig. hi). Corresponding values of X a are 
also shown in the figure. For A a > we used the system size 
L — 1000, and up to 10 5 independent runs were made for each 
curve. The dotted lines have the slope corresponding to the 
exponent S of the dynamical and the directed percolation. 

En Sec. [n] we already noticed the similarity of the 
present model with a forest-fire model with immuniza- 



tion. Results presented in this section provide farther 
arguments supporting such an analogy. Indeed, dynam- 
ical exponents measured by Albano for the forest fire 
model are also close to the dynamical percolation (with- 
out growth of tree) and the directed percolation (with 
growth of tree) JTgJ . 



C. Inhomogeneous absorbing states 



Static simulations suggest that the model becomes crit- 
ical on the line Af, > A^ and A a = 0. Moreover, let us 
notice that for X a = there arc infinitely many absorb- 
ing states: indeed, any configuration without predators 
is an absorbing state. It is well known that dynamical 
Monte Carlo method can be applied also to models with 
infinitely many absorbing states. However, as we will see 
below, applicability of this method to the criticality on 
this line requires serious reconsiderations. 

First, let us recall that the dynamical Monte Carlo 
method for models with infinitely many absorbing states 
usually uses the so-called natural absorbing states, i.e. 
states which are reached by the model's dynamics. Nu- 
merical evidence suggests that for such states the dynam- 
ical critical point coincides with the static one. Moreover, 
dynamical exponents S and rj, measured on such states, 
take universal values. 

Following this prescription, we generated natural ab- 
sorbing states for A a = and Af, > X[ , and then used 
such states to perform dynamical simulations. An initial 
configuration was chosen randomly, with equal probabil- 
ities of a site being empty, occupied by prey, or simulta- 
neously by prey and predator. Fixing A a (= 0) and A;,, 
we then allowed the system to evolve until an absorb- 
ing configuration was reached (i.e. all predators die out). 
Our results, presented in Fig. g, show, however, that the 
spreading of activity is not critical (i.e. power-law), but 
rather exponential. But we can argue that this is not 
surprising. Indeed, a random initial configuration (with 
probability of prey equal 1/3) is below the percolation 
threshold with respect to the clusters of prey, and con- 
tains only finite clusters of them [^9| . On such clusters 
activity has certainly finite duration (for A a = prey do 
not grow) and the exponential decay of P(t), as it can be 
seen in Fig. [| is an expected feature. 

Clearly, the lack of criticality in P(t) is due to the 
finiteness of prey clusters in the natural absorbing states. 
In principle, we can cure this effect starting from ran- 
dom initial configurations but containing larger fraction 
of prey. For sufficiently large concentrations the sys- 
tem will be above the percolation threshold and activity 
will be able to spread infinitely. Probably, for a certain 
concentration of prey we can tune the system to have 
power-law decay for Pit). Such a procedure, however, is 
somehow artificial, and criticality of spreading will not 
be related with the static criticality of the system. 
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FIG. 5. The survival probability P(t) as a function of t ob- 
tained for A a = and for A& = 6 (solid line) and 8 (dotted 
line). Absorbing states were obtained using continuous cool- 
ing with the cooling rates (from the top) r — 0.0001, 0.0001, 
0.001, 0.01, and oo. For the slowest cooling we used the sys- 
tem size L — 1000. For each cooling rate we generated 10 3 
absorbing states and for each absorbing state we generated 
from 10 2 to 10 5 independent runs. 

The question arises here, whether it is possible to gen- 
erate absorbing states on which spreading would be crit- 
ical, and where this criticality would be generated "more 
naturally" ? Hopefully, criticality of spreading on absorb- 
ing states obtained during such a procedure should be re- 
lated with the steady-state criticality of the model. In the 
following we suggest a procedure which imitate the quasi- 
static approach to the critical point on the line A a = 0. 
In our approach we gradually reduce A a according to the 
formula: 

X a = A°exp(-rt), (9) 

where A° = 1 and r is the 'cooling' rate. [We expect that 
the detailed time dependence in Eq. (|J) is not relevant 
as long as it is a slow process] . We terminate the cooling 
when an absorbing state is reached. When the cooling is 
slow, the system have enough time to build large clusters 
of prey. Our simulations for Ab = 6 suggest (see Fig. ||) 
that in the limit r — > 0, such absorbing states are crit- 
ical, with S — 0.59(10) (along trajectories of type 5 in 
Fig. Q). Measuring the number of active sites N(t), and 
using Eq. (jjj) we estimate rj = 0.34(10) (see Fig. |). The 
departure of the curves from a straight line observed for 
large values of time is related to the finiteness of the cool- 
ing rate. Moreover, we measured the averaged squared 
distance R 2 (t), and using Eq. (j|) we obtain z = 2.0(1) 
(see Fig. |?]). Actually, we expect that the correct value 
of this exponent is z = 2. Indeed, in Eq. (j^) one makes 
the average only over surviving runs, thus the long-time 
contributions to R 2 (t) come from the rare events, when 
the activity happened to be placed on a large island of 
prey. Since Af, > A^, the activity on such islands spreads 
in a deterministic way (annular growth), which leads to 
z = 2. 



FIG. 6. The number of active sites N(t) as a function of 
t obtained for X a = 0, At = 6, and 8. Absorbing states 
were obtained using continuous cooling with the cooling rate 
r — 0.0001. Straight dotted lines have slopes corresponding 
to rj = 0.67 and r\ = 0.35. 
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FIG. 7. The squared distanced of active sites R 2 (t) as a 
function of t obtained for A a = 0, Xt — 6 and 8. Absorb- 
ing states were obtained using continuous cooling with the 
cooling rate r = 0.0001. The straight dotted line has slope 
corresponding to z = 2. 



The same procedure for Ab = 8 yields: 5 = 0.35(10), 
r) = 0.67(10) and z — 2.0(1). (Relatively large errors of 
estimations of critical exponents are due to several, and 
difficult to estimate factors, as finite time t, finite cooling 
rate r, and statistical fluctuations.) Such results confirm 
that z = 2, and suggest that exponents S and 77 might 
change continuously with Ab- Non- universality of these 
exponents is a well-known property for some other mod- 
els belonging to the directed percolation universality class 
|p0f . Note, however, that the situation is different in our 
case, because the non-universality is related to the value 
of Ab, rather than to the choice of the initial state. Such 
control parameter dependence has already been observed 
in other models |21f|. Note that non-universal behavior 
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is not present along the DP line since the corresponding 
absorbing state is unique. Let us finally note that 8 + 77 
seems to be close to unity, which is an exact mean-field 
result (5mf = 1,Vmf = 0). This is the only dynamical 
trace of the mean-field nature of the transition observed 
in the steady state. Let us emphasize, however, that 
the criticality of spreading appears only if we prepare 
the absorbing states using the method which mimics the 
quasi-steady-state evolution of the model. 

V. CONCLUSIONS 

The detailed investigation of the critical properties of 
a simple prey-predator model introduced in Ref. flC[ | 
showed the presence of three different types of non- 
equilibrium phase transitions between active and absorb- 
ing states. First, the existence of a usual DP-like transi- 
tion line was confirmed at AJ(A a ) for A a > 0. Second, 
a mean-field-like transition was observed for A a — » 0, 
Ab > X[ . The mean-field character of this transition 
can be explained in terms of oscillations present in the 
model. Indeed, as described in Ref. Jl0| , when X a — > 0, 
the system is subject to large density oscillations. These 
oscillations are generating an important local mixing of 
the possible states, leading to a mean-field-like behavior. 
The criticality along the A a = line was confirmed with 
the dynamical approach using specially prepared inho- 
mogeneous initial states. Some dynamical trace of the 
mean-field nature of this transition was also observed. 
Third, at the bicritical point T, where the two different 
critical lines meet, we found a dynamical percolation type 
transition moving along the A a = line, while approach- 
ing the T point from finite A a 's we observed a new type 
of critical behavior. 

The measured exponents corresponding to the A a = 
line are summarized in Table |[ The best numerical es- 
timates for the critical exponents of the two-dimensional 
dynamical percolation are given for comparison: <5 = 
0.092, r\ = 0.586, and z = 1.771 g§. 
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A b = 4.67 


5.0 


6.0 8.0 


/3 2 


1.33(4) 


1.01(1) 


0.96(4) 
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-0.65(10) 


-0.1(1) 


-0.05(5) 
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0.095(5) 




0.59(10) 0.35(10) 
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0.60(5) 
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1.72(4) 
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TABLE I. Critical exponents around A a = 
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